1. Cel

Celem przedmiotowego ćwiczenia jest zapoznanie się studentów z pracą z danymi pochodzącymi z bazy GIOś oraz z danymi meteorologicznymi. Aby prawidłowo wykonac zlecone nam zadanie musielismy zapoznac sie z formatem danych, poznać narzędzia które ułatwią nam pracę oraz ich przetwarzanie. Musieliśmy dokonać prognozy wybranego przez nasz zanieczyszczenia powietrza na dany rok, oraz dokonać analizy w celu określenia największej korealcji danych aby z łatwością móc odczytać panujące trendy w wybranej stacji obserwacyjnej. Ostatnim aspektem naszej pracy, było napisanie sprawozdania technicznego z wykonanengo ćwiczenia przy pomocy RMarkdown, w celu udokumentowania naszej pracy w języku R.

2. R (język programowania)

R - interpretowany jezyk programowania oraz srodowisko do obliczen statystycznych. Stosowany jest w analizie szeroko rozumianych danych środowiskowych i przestrzennych oraz ich wizualizacji. Podobny jest do języka i środowiska S stworzonego w Bell Laboratories przez Johna Chambersa i jego współpracowników. R jako implementacja języka S została stworzona przez Roberta Gentlemana i Rossa Ihakę na uniwersytecie w Auckland. Nadaje się on świetnie do interaktywnej pracy z danymi, ponieważ połączono w nim wybrane cechy języków funkcyjnych oraz obiektowych.

Źródlo

3. R Markdown

R Markdown - jest formatem pliku stworzonym do sporzadzania dynamicznych dokumentow z wykorzystaniem R. Plik typu Markdown jest pisany w specyficzny dla siebie sposob, który zaklada bardzo latwa edycje tekstu oraz implementowanie w nim fragmentow kodu (chunki zawierajace kod z poleceniami w jezyku R). R Markdown jest bardzo wygodna metoda formatowania plikow HTML, PDF i dokumentow MS Word.

Źródlo

4. Baza GIOŚ

Baza GIOŚ - (Główny Inspektorat Ochrony Środowiska) - to zbiór publicznie dostępnych danych zebranych przez organ jakim jest GIOŚ.

Jednym z najistotniejszych zadań realizowanych przez Inspekcję Ochrony Środowiska jest prowadzenie badań i ocen stanu środowiska, w tym monitoringu jakości powietrza. Zadanie to jest wykonywane w ramach Państwowego Monitoringu Środowiska (PMŚ), którego program jest opracowywany przez Głównego Inspektora Ochrony Środowiska i zatwierdzany przez Ministra Środowiska. W oparciu o krajowy program PMŚ opracowywane są wojewódzkie programy PMŚ zatwierdzane przez Głównego Inspektora Ochrony Środowiska.

Monitoring jakości powietrza obejmuje zadania związane z badaniem i oceną stanu zanieczyszczenia powietrza, w tym pomiary i oceny jakości powietrza w strefach, monitoring tła miejskiego pod kątem WWA, pomiary stanu zanieczyszczenia powietrza pyłem PM2,5 dla potrzeb monitorowania procesu osiągania krajowego celu redukcji narażenia, pomiary stanu zanieczyszczenia powietrza metalami ciężkimi i WWA oraz rtęcią w stanie gazowym na stacjach monitoringu tła regionalnego, pomiary składu chemicznego pyłu PM2,5, monitoring prekursorów ozonu; programy badawcze dotyczące zjawisk globalnych i kontynentalnych wynikające z podpisanych przez Polskę konwencji ekologicznych.

Ok. 90% pomiarów jakości powietrza wykonywanych w ramach PMŚ oraz roczne i pięcioletnie oceny jakości powietrza w strefach są wykonywane przez GIOŚ i regionalne wydziały. Na zlecenie Głównego Inspektoratu Ochrony Środowiska są realizowane krajowe programy monitoringu jakości powietrza, GIOŚ jednocześnie nadzoruje i koordynuje wykonywanie programu badań i ocen jakości powietrza określonego w krajowym i wojewódzkich programach Państwowego Monitoringu Środowiska.

Źródlo

5. Las Losowy (Random Forest)

Las Losowy (ang. Random Forest) - jedna z metod ensamblingu uczenia maszynowego. Jest to rozszerzenie metody drzew decyzyjnych, która polega na budowaniu modelu w oparciu o zadawanie i dzielenie zbioru według warunków nakładanych na dane zmienne. I chociaż lasy losowe nadają się lepiej do prognozowania danych dyskretnych bądź nominalnych (metoda klasyfikacyjna), to również w naszym przypadku powinna przynieść satysfakcjonujące efekty.

[Źródlo] Materiał Zajęciowy

6. Wybrana Stacja

Aby rozpocząc pracę z naszymi danymi, musieliśmy je pozyskać, w tym celu skorzystaliśy z środowiska dołączonego przez Pana Doktora, jako materiał zajęciowy “projekt_2.RData”. Aby zaimportować dane skorzystalićmy z poniższego polecenia:

load(file = "projekt_2.RData")

Naszą pracę z danymi postanowiliśmy rozpocząć od usnięcia polskich znaków w kolumnach z których w dalszczej części projektu będziemy korzystać:

colnames(gios_inv)[c(7,13)] <- c("Data.zamkniecia", "Miejscowosc")
colnames(gios_inv)[15:16] <- c("Szerokosc", "Dlugosc")

W celu wykonania naszego ćwiczenai postanowiliśmy skorzystać z stacji znajdującej się północnej cześci Częstochowy, w dzielnicy “Tysiąclecia” na ulicy. Krzysztofa Kamila Baczyńskiego. W tym celu skorzystaliśmy z stacji SlCzestoBacz.

gios_inv %>% filter(Kod.stacji == "SlCzestoBacz") -> czesto_lok

czesto_lok_info <- paste(paste(
  czesto_lok$Kod.stacji,
  paste("Miejscowość:", czesto_lok$Miejscowosc),
  paste("Data uruchomienia:", czesto_lok$Data.uruchomienia),
  paste("Data zamknięcia:", czesto_lok$Data.zamkniecia),
  paste("Typ stacji:", czesto_lok$Typ.stacji),
  paste("Typ obszaru:", czesto_lok$Typ.obszaru),
  paste("Współrzędne:", czesto_lok$Dlugosc,"E, ", czesto_lok$Szerokosc,"N"),
  sep = "<br/>"
))


leaflet() %>%
  addTiles() %>%
  addMarkers(data=czesto_lok,
             lng= ~ Dlugosc,
             lat= ~ Szerokosc,
             popup = czesto_lok_info)

Lokalizację geograficzną naszej stacji jakości powietrza sprawdziliśmy nieco wcześniej aby być pewnymi że w jej poblizu znajduje się stacja meteorologiczna z której będziemy mogli skorzystać dla prawidłowej analizy. Jak widzimy na poniższej mapie, odległość dzieląca nasze stacje nie przekracza 4km, dzięki czemu dane będą bardzo obiektywne.

noaa_isd <- getMeta(end.year="current", lon=19.130111, lat=50.836389, returnMap=T)

noaa_isd

Po znalezieniu interesujacej nas stacji, postnaowiliśy zaimportować interesujące nas dane do zmiennej, tak aby w dalszej częsci naszej pracy nad projektem, z łatwością móc z nich korzystać, bez wracanie do tego momentu przygotowywania danych.

czestochowa_met <- importNOAA(code="125500-99999", year=2000:2020)

Postanowiliśmy dokonać analizy ciagłości obserwacji, dla wybranej przez nas stacji, aby tego dokonać skorzystaliśmy z poniższego fragmentu kodu:

# Dodanie informacji o roku do ramki danych

czestochowa_met$years <- format(czestochowa_met$date, "%Y")

# Policzenie ilości pomiarów dla danego roku dla stacji Okęcie

czestochowa_met %>%
  group_by(years) %>%
  summarise(liczba_pomiarow = n() - sum(is.na(wd))) -> czestochowa_met_summary

czestochowa_met_summary <- czestochowa_met_summary %>% 
  as_tibble() %>% 
  mutate(years = as.integer(years))

czy_pomiary_czestochowa_met <- data.frame(ymin = -Inf,
                                 ymax = Inf,
                                 xmin = c(1999 , 2021),
                                 xmax = c(1999.5 , 2021.5))

# Utworzenie wykresu

ggplot(czestochowa_met_summary, aes(x=years, y=liczba_pomiarow, color=liczba_pomiarow)) +
  geom_point(size=3, show.legend = FALSE) + 
  scale_y_continuous(name = 'Liczba rekordów pomiarowych',
                     limits = c(0, 9000),
                     expand = c(0, 0)) +
  scale_x_continuous(name = 'Rok',
                     limits = c(1998, 2022),
                     expand = c(0, 0),
                     breaks = seq(2000, 2020, by = 1)) +
  geom_rect(data = czy_pomiary_czestochowa_met, 
            aes(ymin = ymin, ymax = ymax,
                xmin = xmin, xmax = xmax), fill = 'red',
            alpha = 0.15, inherit.aes = FALSE) +
  theme_minimal() +
  theme(plot.margin = unit(c(.5,.5,.5,.5), "cm")) +
  labs(title = 'Liczba pomiarów dla stacji pomiarowej CZESTOCHOWA (CZESTOCHOWA) w ISD NOAA', 
       subtitle = paste0('Stan na ',format(Sys.time(),'%d %B, %Y')))

Stacja Częstochowa działa od 1940 roku, w badanym przez nas okresie czasu - to jest od 2000 - 2020 roku, stacja nie posiadała przerwy w dostarczaniu danych meteorologicznych. dodatkowo z powyższego wykresu możemy zauważyć, że najlepszą dokładnością danych charakteryzuje się okres od 2017-2019, gdyż to właśnie w tym okresie wykonywanych było najwięcej pomiarów, oraz zebrane dane miały charakter ciągły.

7. Analiza danych

7.1. Selekcja danych

Analizę naszych danych rozpoczeliśmy od wyselkecjonowania interesujących nas danych dotyczących zanieczyszczenia pyłami PM10 z wybranej przez nas stacji. Pozwoli nam to znacząco ograniczyć czas wykonywanych poleceń, oraz pozwoli skupić się na dokładniejszej analizie wybranego zagadnienia.

czestochowa_PM10 <- PM10_1h %>%
  filter(kod=="SlCzestCzes_baczy" | kod=="SlCzestoBacz") %>% 
  select(date, obs)

Aby nasze dane prawidłowo z sobą współgrały musieliśy wykonać szereg poleceń mających na celu dodanie godziny z uwzględnieniem stefy czasowej, dokonać selekcji wybranych danych a następnie połaczyć je w jeden plik z danymi. W tym celu skrozystaliśmy z poleceń :

## Musimy dodać godzinę w danych meteo, aby uwzględnić różne strefy czasowe danych
czestochowa_met$date <- czestochowa_met$date + 3600

## selekcja wybranych zmiennych  

czestochowa_met <- czestochowa_met[c(3, 7:14)]

## Łączenie danych meteorologicznych z danymi jakości powietrza
dane <- inner_join(czestochowa_met, czestochowa_PM10, by = "date")

Nasze dane dla lepszego zobrazowania postanowiliśmy przedstawić w formie tabeli danych aby każdy mógł odczytać interesujące go dane

datatable(dane)

W rozpatrywanym okresie występuje wysoka kompletnosć danych. W dalszej czesci naszej analizy posluzymy sie wynikami w formie graficznej, gdyz przy takiej ilosci danych, taka forma bedzie bardziej czytelna.

7.2. Braki danych i selekcja

Wykorzystując to metodę niezbędne jest posiadanie kompletnych wierszy obserwacji, a także tylko tych zmiennych, które wniosą informację dodaną do naszego modelu. Aby upewnić się w którym roku posiadanmy najlepszą kompletność danych, skorzystaliśmy z poniższego polcenia, jednocześnie wybierając te które charakteryzowały się najlepszą kompletnością.

# Sprawdzenie braków danych i poszczególnych kolumn
summary(dane)
##       date                           ws              wd           air_temp     
##  Min.   :2006-01-01 00:00:00   Min.   : 0.00   Min.   :  0.0   Min.   :-25.70  
##  1st Qu.:2009-04-01 20:00:00   1st Qu.: 2.00   1st Qu.:140.0   1st Qu.:  2.20  
##  Median :2012-07-01 17:00:00   Median : 2.00   Median :210.0   Median :  9.00  
##  Mean   :2012-07-01 20:36:55   Mean   : 2.48   Mean   :201.6   Mean   :  9.06  
##  3rd Qu.:2015-10-02 00:00:00   3rd Qu.: 3.00   3rd Qu.:280.0   3rd Qu.: 16.00  
##  Max.   :2018-12-31 23:00:00   Max.   :13.00   Max.   :360.0   Max.   : 36.50  
##                                NA's   :42861   NA's   :43668   NA's   :42216   
##    atmos_pres      visibility      dew_point            RH        
##  Min.   : 976    Min.   :    0   Min.   :-28.60   Min.   : 13.03  
##  1st Qu.:1012    1st Qu.: 5000   1st Qu.: -0.50   1st Qu.: 65.10  
##  Median :1017    Median :10000   Median :  5.00   Median : 81.98  
##  Mean   :1017    Mean   :12995   Mean   :  4.61   Mean   : 77.00  
##  3rd Qu.:1022    3rd Qu.:20000   3rd Qu.: 10.50   3rd Qu.: 92.38  
##  Max.   :1051    Max.   :50000   Max.   : 21.90   Max.   :100.00  
##  NA's   :42250   NA's   :42316   NA's   :42250    NA's   :42250   
##     ceil_hgt          obs        
##  Min.   :   15   Min.   :  0.00  
##  1st Qu.:  750   1st Qu.: 15.36  
##  Median : 1200   Median : 25.43  
##  Mean   : 8841   Mean   : 35.81  
##  3rd Qu.:22000   3rd Qu.: 43.00  
##  Max.   :22000   Max.   :705.91  
##  NA's   :92012   NA's   :13122
# Selekcja tylko kompletnych wierszy
dane_rf <- dane_rf[complete.cases(dane_rf),]

W celu przeprowadzenia dokładniejszej analizy postanowiliśmy przeprowadzić analizę dla więcej niż jednego modelu. W tym celu postanowiliśmy utworzyć nowe zmienne naszych danych różniące się parametrami, aby otrzymać różny punkt odniesienia, na którym będziemy mogli bazować. Zdecydowaliśmy się na utworzenie 5 zmiennych gdyż taka ilość nie powodowała błędów wykonywanego kodu, oraz czas pracy wykonywanej przez system operacyjny był stosunkowo szybki.

dane_rf %>% mutate(years = year(date)) -> dane_rf

#dodanie kolumn
dane_rf$unix_date <- as.numeric(as.POSIXct(dane_rf$date))

dane_rf %>% mutate(years = year(date)) -> dane_rf

dane1 <- dane_rf %>%
  mutate(week = week(date))

dane2 <- dane1 %>% 
  mutate(hour = hour(date))

dane3 <- dane2 %>% 
  mutate(weekday = wday(date))

dane4 <-dane3 %>% 
  mutate(month = month(date))

Wykreślmy wykres, który poznaliśmy w czasie wykonywania Projektu numer jeden, który ma na celu zobrazowanie ilości rekordów w danym roku:

# Policzenie ilości pomiarów w poszczególnych latach
dane_rf$years <- format(dane_rf$date, "%Y")

dane_rf %>%
  group_by(years) %>%
  summarise(liczba_pomiarow = n()) -> dane_rf_summary

ggplot(dane_rf_summary, aes(x=years, y=liczba_pomiarow, color=liczba_pomiarow)) +
  geom_point(size=3, show.legend = FALSE) +
  theme_minimal() +
  coord_flip()

Zgodnie z przykładowym projektem II, dołączonym jako materiał zajęciowy, dzielimy wszystkie zbiory na dwa podzbiory. Jeden z nich, tzw. “zbiór treningowy”, będzie obejmować dane, na których nauczymy model przewidywać stężenia PM10 na podstawie danych meteorologicznych. Nie wybierajmy zbyt wiele rekordów, bo zwiększy to wykładniczo czas obliczeniowy. Zdecydujmy się na 2 - 3lata najbardziej kompletnych danych. My do zbioru treningowego postanowiliśmy skorzystać z lat 2016 - 2018. Do zbioru, testowego, czyli takiego, na którym będziemy weryfikować nasz model, wybieramy rok 2018. Właśnie dla tego przedziału skorzystaliśmy dla każdej z naszych zmiennych.

# Filtrujemy dane do zbioru treningowego i testowego

dane_rf %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane_rf_train

dane_rf %>%
  filter(years==2018) %>%
  select(-years) -> dane_rf_test


dane1 %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane1_rf_train

dane1 %>%
  filter(years==2018) %>%
  select(-years) -> dane1_rf_test


dane2 %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane2_rf_train

dane2 %>%
  filter(years==2018) %>%
  select(-years) -> dane2_rf_test


dane3 %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane3_rf_train

dane3 %>%
  filter(years==2018) %>%
  select(-years) -> dane3_rf_test


dane4 %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane4_rf_train

dane4 %>%
  filter(years==2018) %>%
  select(-years) -> dane4_rf_test

Aby w przyszłosći nie pojawił się problem spójności danych, postanowiliśmy zapisać wyselekcjonowane dane jako nowa zmienna. W dalszej analizie bedziemy bazować na, kopii tak aby nie uszkodzić pliku oryginalnego.

# Zapisanie kopii obiektu, tak aby nie zmieniać nazw kolumn w oryginalnym zestawie danych
dane_rfc <- dane_rf_train
dane1_rfc <- dane1_rf_train
dane2_rfc <- dane2_rf_train
dane3_rfc <- dane3_rf_train
dane4_rfc <- dane4_rf_train

7.3. Korelacje między zmiennymi

Aby sprawdzić korelacje liniowe między naszymi zmiennymi posłużymy się korelacją Pearsona - tą samą z której krozystaliśmy w Projekcie I.

Macierz korelacji liniowej Pearsona (macierz współczynników określających poziom zależności liniowej między zmiennymi losowymi), pozwala wstępnie sprawdzić występowanie zależności między analizowanymi zmiennymi. Mówiąc prosto, współczynnik ten może przyjmować wartości od -1 do 1. Im jest mniejszy, tym silniejsza zależność ujemna między parametrami, a im jest bliższy 1, tym silniejsza zależność dodatnia. Dodatkowo możemy sprawdzić, czy otrzymane wartości są istotne statystycznie (przyjmiemy poziom istotności równy 0.05).

Aby mówić w ogóle o wykorzystaniu modelu, musimy mieć solidne podstawy, aby sądzić, że związek między zmiennymi jest silny i jedną z nich (PM10 - zmienna objaśniana) da się przedstawić jako funkcja pozostałych zmiennych (zmienne objaśniające).

# Definiowanie indywidualnej palety kolorów w zależności od współczynnika korelacji (kodowanie hex)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

# Policzenie odpowiednich współczynników za pomocą funkcji 'corr'
corr <- rcorr(as.matrix(dane_rfc))

# Zdefiniowanie wektora ze współczynnikami korelacji
corr_r <- corr$r

# zdefiniowanie wektora z wartościami p-value
corr_p <- corr$P

# Wykreślenie macierzy korelacji
corrplot(corr_r, method = "color", col = col(200),  
         type = "upper", order = "hclust", addCoef.col = "black", 
         diag = FALSE,
         tl.col = "black", 
         tl.srt = 45, 
         p.mat = corr_p, 
         sig.level = 0.05)

Jak widzimy, korelacje liniowe Pearsona między stężeniem pyłu PM10 a innymi zmiennymi są jedynie na umiarkowanym poziomie, ale wszystkie są statystycznie istotne. Ta wiedza jest dla nas wystarczająca, aby być pewnym że na naszych danych możemy zbudować model, jednocześnie nie musimy przejmować się, że jakaś korealacja okaże się nieistotna statystycznie - model prawdopodobnie sam ją odrzuci na etapie uczenia.

7.4. Wykonanie predykcji

Do uruchomienia lasu losowego wykorzystano pakiet caret, który jest bardzo bogatym pakietem do modelowania statystycznego w R. Aby oszczędzić nieco czasu zgodnie z propozycją zawartą w konspekcie uruchomimy także obliczenie równoległe, które wykona się wykorzystując dostępne rdzenie procesora. Skorzystaliśmy z takiej metody gdyż utworzenie 5 modeli, było bardzo czasochłonne na naszym sprzęcie, a ta metoda pozwoliła odrobinę oszczędzić nam czasu, oraz ograniczyłą zużycie zasobów naszego sprzętu.

# Ustawiamy metodę tzw. "kroswalidacji", czyli samo-testowania się modelu w trakcie jego budowy.
cross.walid <- trainControl(method = "cv",
                            number = 5,
                            allowParallel = T,
                            returnResamp = 'all')

# Zadajemy, aby w modelu były obecne wszystkie zmienne oprócz objaśnianej (PM10)
tunegrid <- expand.grid(.mtry = c(1:10))
tunegrid1 <- expand.grid(.mtry = c(1:11))
tunegrid2 <- expand.grid(.mtry = c(1:12))
tunegrid3 <- expand.grid(.mtry = c(1:13))
tunegrid4 <- expand.grid(.mtry = c(1:14))
tunegrid5 <- expand.grid(.mtry = c(1:15))
# Uruchamiamy klaster obliczeniowy z wszystkich dostępnych rdzeni - 1 (zwyczajowo zostawia się jeden, aby system swobodnie działał)
cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)

# Budujemy (trenujemy) model lasu losowego
model_rf <- train(obs ~ ., 
                  data = dane_rf_train, 
                  method = "rf", 
                  replace = T,
                  ntree = 200,
                  metric = 'RMSE',
                  tuneGrid = tunegrid, 
                  trControl = cross.walid)

model1_rf <- train(obs ~ ., 
                  data = dane1_rf_train, 
                  method = "rf", 
                  replace = T,
                  ntree = 200,
                  metric = 'RMSE',
                  tuneGrid = tunegrid1, 
                  trControl = cross.walid)

model2_rf <- train(obs ~ ., 
                   data = dane2_rf_train, 
                   method = "rf", 
                   replace = T,
                   ntree = 200,
                   metric = 'RMSE',
                   tuneGrid = tunegrid2, 
                   trControl = cross.walid)

model3_rf <- train(obs ~ ., 
                   data = dane3_rf_train, 
                   method = "rf", 
                   replace = T,
                   ntree = 200,
                   metric = 'RMSE',
                   tuneGrid = tunegrid3, 
                   trControl = cross.walid)

model4_rf <- train(obs ~ ., 
                   data = dane4_rf_train, 
                   method = "rf", 
                   replace = T,
                   ntree = 200,
                   metric = 'RMSE',
                   tuneGrid = tunegrid4, 
                   trControl = cross.walid)




# Zatrzymujemy klaster
stopCluster(cluster)
registerDoSEQ()

Z poniższych wykresów możemy z łatwością zobaczyć zmiany w testowaniu się modelu.

#Model 9 zmiennych
plot(model_rf)

#Model 10 zmiennych
plot(model1_rf)

#Model 11 zmiennych
plot(model2_rf)

#Model 12 zmiennych
plot(model3_rf)

#Model 13 zmiennych
plot(model4_rf)

Jednak aby dokładnie odczytać wartość mtry, RMSE, Rsquared oraz MAE, skorzystamy z poniższych wyników .

#Model 9 zmiennych
model_rf
## Random Forest 
## 
## 11184 samples
##     8 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8947, 8947, 8947, 8947, 8948 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##    1    22.14724  0.8040095  12.94263
##    2    20.82950  0.8196639  12.17065
##    3    20.58272  0.8228228  12.03239
##    4    20.50626  0.8232116  11.96935
##    5    20.49883  0.8227717  11.95867
##    6    20.45653  0.8231178  11.94832
##    7    20.59064  0.8205505  12.02360
##    8    20.75866  0.8173084  12.09947
##    9    20.70904  0.8182752  12.08131
##   10    20.70458  0.8183340  12.07470
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 6.
#Model 10 zmiennych
model1_rf
## Random Forest 
## 
## 11184 samples
##    10 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8948, 8946, 8948, 8947, 8947 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##    1    19.53636  0.8493156  11.33375
##    2    18.16269  0.8648520  10.50227
##    3    18.06996  0.8647392  10.41039
##    4    18.13147  0.8631444  10.43246
##    5    18.17947  0.8616916  10.44163
##    6    18.22309  0.8609007  10.45905
##    7    18.23165  0.8605247  10.51912
##    8    18.37433  0.8582054  10.55277
##    9    18.35812  0.8581134  10.60079
##   10    18.58768  0.8544083  10.69067
##   11    18.58530  0.8542627  10.69724
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 3.
#Model 11 zmiennych
model2_rf
## Random Forest 
## 
## 11184 samples
##    11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8947, 8948, 8948, 8946, 8947 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##    1    19.64163  0.8515474  11.46933
##    2    17.82246  0.8707086  10.30115
##    3    17.60785  0.8722647  10.20005
##    4    17.59252  0.8721022  10.19026
##    5    17.64774  0.8705475  10.22497
##    6    17.70192  0.8695952  10.21369
##    7    17.66486  0.8699411  10.23112
##    8    17.84399  0.8669006  10.30405
##    9    17.77658  0.8679943  10.31162
##   10    17.86031  0.8660553  10.37553
##   11    18.05088  0.8630838  10.44521
##   12    17.96967  0.8642814  10.41103
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
#Model 12 zmiennych
model3_rf
## Random Forest 
## 
## 11184 samples
##    12 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8948, 8947, 8946, 8947, 8948 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##    1    19.32859  0.8604909  11.25531
##    2    17.52587  0.8771690  10.18301
##    3    17.30842  0.8778948  10.05658
##    4    17.19356  0.8783876  10.00094
##    5    17.31842  0.8760826  10.04116
##    6    17.32375  0.8754202  10.06966
##    7    17.32627  0.8748956  10.09337
##    8    17.50574  0.8720442  10.14486
##    9    17.53227  0.8714417  10.15545
##   10    17.52422  0.8713014  10.19491
##   11    17.62599  0.8694167  10.24135
##   12    17.69618  0.8681833  10.28774
##   13    17.65027  0.8685643  10.28715
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.
#Model 13 zmiennych
model4_rf
## Random Forest 
## 
## 11184 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8947, 8948, 8948, 8947, 8946 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE      
##    1    19.81016  0.8546864  11.395945
##    2    17.61938  0.8761540  10.080288
##    3    17.44280  0.8760401   9.939129
##    4    17.35046  0.8767173   9.907598
##    5    17.40871  0.8753637   9.948455
##    6    17.52787  0.8730406  10.006474
##    7    17.51726  0.8731641  10.017711
##    8    17.58470  0.8716453  10.070079
##    9    17.66427  0.8703834  10.105871
##   10    17.65625  0.8705696  10.129317
##   11    17.79102  0.8683707  10.196462
##   12    17.79616  0.8682362  10.220303
##   13    17.95914  0.8654848  10.306924
##   14    17.94380  0.8657578  10.313816
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.

7.5. Testowanie modelu

Wykorzystując przed chwilą storzone modele możemy sprawdzić jak zachowują się one na danych, które nie były wykorzystane do jego uczenia. Pobierane zostaną dane meteorologiczne dla roku 2018 dla tej samej stacji w Częstochowie, a następnie wykorzystany zostanie utworzony model do oszacowania stężenia pyłu PM10 w tym roku.

dane_rf_test$mod <- predict(model_rf, dane_rf_test)
dane1_rf_test$mod <- predict(model1_rf, dane1_rf_test)
dane2_rf_test$mod <- predict(model2_rf, dane2_rf_test)
dane3_rf_test$mod <- predict(model3_rf, dane3_rf_test)
dane4_rf_test$mod <- predict(model4_rf, dane4_rf_test)

Żeby sprawdzić, jak dobrze nasz model odwzorowuje rzeczywistość musimy mieć jakąś wartość odniesienia. Sporządźmy dwa wykresy obrazujące nam zależność wartości stężenia PM10 zaobserwowanego na stacji jakości powietrza od zamodelowanej jego wartości. Pierwszy z nich to wykres rozrzutu.

# Wykres rozrzutu

# 9 zmiennych
ggplot(data=dane_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
  theme_minimal()

# 10 zmiennych
ggplot(data=dane1_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
  theme_minimal()

# 11 zmiennych
ggplot(data=dane2_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
  theme_minimal()

# 12 zmiennych
ggplot(data=dane3_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
  theme_minimal()

# 13 zmiennych
ggplot(data=dane4_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Czestochowie w roku 2018") +
  theme_minimal()

Wykres rozrzutu pokazuje jak mocno wyniki prognoz odchylone są od pomiarów stężeń PM10. Szara przerywana linia obrazuje model idealny (w sytuacji, gdy prognoza równa jest pomiarowi, czyli x=y), a niebieska – linię trendu uzyskaną z naszych danych, o równaniu wyświetlonym na wykresie.

Według naszych danych, prognoza jest zaniżona w stosunku do pomiarów stężenia PM10 na stacji w Częstochowie, co jest szczególnie widoczne w przypadku najwyższych stężeń. Może to wynikać z faktu, że warunki meteorologiczne nie są jedynymi czynnikami wpływającymi na stężenie zanieczyszczenia na danym obszarze. Pominęliśmy całkowicie kwestię emisji, która, w przypadku sytuacji wyjątkowych i chwilowych, może drastycznie wpłynąć na mierzone stężenia zanieczyszczeń.

Uznaliśmy że porównanie tych wykresów kiedy znajdują się jeden pod drugim, jest bardzo nie praktyczne. W tym celu postanowiliśy przedstawić je bardziej czytelnie. W tym celu postanowiliśy opisać każdy z wykresóW zgodnie z jego modelem, a następnie zrobić jeden obraz zbiorczy, w tym celu skrozystalisy z poniższych poleceń:

bind_rows(dane_rf_test %>% 
             mutate(typ = 'Model A [+unix_date]'),
           dane1_rf_test %>% 
             mutate(typ = "Model B [+week]"),
           dane2_rf_test %>% 
             mutate(typ = "Model C [+hour]"),
           dane3_rf_test %>% 
             mutate(typ = "Model D [+weekday]"),
           dane4_rf_test %>% 
             mutate(typ = "Model E [+month]")) -> new_data
        


ggplot(data=new_data, aes(x=obs, y=mod))+
  geom_point(alpha=0.6, color="purple") +
  geom_smooth(method="lm", formula = y~x-1) +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nd stacji w Czestochowie w roku 2018") +
  theme_minimal() +
  facet_wrap(~typ)

Dzięki takiej formie prezentacji możemy, z łatwością porównać utworzone modele, odczytać interesujace nas parametry, oraz określić który model charakteryzuje się najlepszą dokładnoscią.

Wykresy dodatkowo dostarcza informacji o równaniu regresji to znaczy o równaniu linii dopasowanej do danych przedstawionych na wykresie. Równanie dla najlepszego wyniku naszego doświadczenia przedstawia się następująco: y = 0, 841x.

Dla danych powyższych wykresu wyliczono też wartość współczynnika determinacji oznaczonego jako R2. Współczynnik ten w przypadku naszych wykresów jest równy = 0,81 (w kilku przypadkach). Zgodnie z przedstawioną poniżej tabelą można wyciagnąć wnioski że rozrzut naszej prognozy jest dobrze dopasowany do pomiaru idealnego.

R2
## # A tibble: 5 x 2
##   X1        X2                         
##   <chr>     <chr>                      
## 1 0.0 – 0.5 dopasowanie niezadowalające
## 2 0.5 – 0.6 dopasowanie słabe          
## 3 0.6 – 0.8 dopasowanie zadowalające   
## 4 0.8 – 0.9 dopasowanie dobre          
## 5 0.9 – 1.0 dopasowanie bardzo dobre

Ostantim aspektem naszej pracy jest utworzenie wykresu trendu prognozy dla obserwacji steżeń PM10. Postanowiliśmy utworzyć jeden wykres, który dzięki wysokiej dokładności, powinen być najbardziej dokłądną reprezentacją danych. Skorzystalśmy z zmiennej dane4_rf_test, gdyż to właśnie dla tych danych wyrkes rozrzuty charakteryzował się najlepszą dokładnością.

# Porównanie serii danych w danym fragmencie roku
ggplot(data = dane4_rf_test) +
  geom_path(aes(x = date, y = obs, color = "Obserwacje"), size = 1.5) + 
  geom_path(aes(x = date, y = mod, color = "Model"), size = 1.5, alpha=0.6) +
  scale_x_datetime(limits = ymd_h(c("2018-01-01 00", "2018-01-31 23"))) +
  scale_y_continuous(limits=c(0,250)) +
  scale_color_manual(values=c("Obserwacje"="grey","Model"="red"))+
  xlab("Data") +
  ylab("Stężenie PM10 [µg/m3]") +
  ggtitle("Wykres trendu prognozy od obserwacji dla stężeń PM10 na stacji w Częstochowie w roku 2018") +
  theme_minimal() +
  theme(legend.position="top",legend.title=element_blank())

Wykres trendu prognoz (linia czerwona) od obserwacji stężeń PM10 (linia szara) pozwala na ogólną ocenę zgodności tych dwóch serii danych. Możemy zauważyć, że wniosek wyciągnięty na podstawie poprzedniego wykresu jest potwierdzony – najwyższe wartości stężeń są zaniżane w naszym modelu. Ogólny trend jest zachowany, chociaż również dla najniższych stężeń zdarzają się odchylenia.

Oprócz grafiki wizualizacji oceny dokąłdności modelu zawsze warto przyjrzeć się paramtrom statystycznym oceny dokładności prognoz.

# 9 zmiennych
modStats(mydata = dane_rf_test, 
         mod = "mod", 
         obs = "obs", 
         type = "season") %>% 
  knitr::kable(digits = 2)
season n FAC2 MB MGE NMB NMGE RMSE r COE IOA
2 spring (MAM) 1413 0.84 -9.37 17.77 -0.18 0.34 28.85 0.79 0.41 0.70
3 summer (JJA) 1185 0.85 -3.62 8.51 -0.15 0.35 12.07 0.42 0.10 0.55
1 autumn (SON) 1276 0.82 -5.69 16.22 -0.13 0.37 24.08 0.56 0.24 0.62
4 winter (DJF) 1945 0.85 0.18 20.02 0.00 0.37 31.52 0.70 0.39 0.69 
# 10 zmiennych
modStats(mydata = dane1_rf_test, 
         mod = "mod", 
         obs = "obs", 
         type = "season") %>% 
  knitr::kable(digits = 2)
season n FAC2 MB MGE NMB NMGE RMSE r COE IOA
2 spring (MAM) 1413 0.88 -5.80 17.28 -0.11 0.33 29.16 0.77 0.42 0.71
3 summer (JJA) 1185 0.81 4.22 9.73 0.17 0.40 12.70 0.37 -0.03 0.49
1 autumn (SON) 1276 0.86 -0.55 15.80 -0.01 0.36 22.25 0.60 0.26 0.63
4 winter (DJF) 1945 0.79 6.60 22.02 0.12 0.40 31.88 0.70 0.33 0.66 
# 11 zmiennych
modStats(mydata = dane2_rf_test, 
         mod = "mod", 
         obs = "obs", 
         type = "season") %>% 
  knitr::kable(digits = 2)
season n FAC2 MB MGE NMB NMGE RMSE r COE IOA
2 spring (MAM) 1413 0.89 -5.97 16.99 -0.12 0.33 28.72 0.78 0.43 0.72
3 summer (JJA) 1185 0.82 3.68 9.49 0.15 0.39 12.49 0.38 0.00 0.50
1 autumn (SON) 1276 0.87 -1.00 15.36 -0.02 0.35 21.96 0.61 0.28 0.64
4 winter (DJF) 1945 0.81 5.96 21.21 0.11 0.39 31.17 0.72 0.35 0.68 
# 12 zmiennych
modStats(mydata = dane3_rf_test, 
         mod = "mod", 
         obs = "obs", 
         type = "season") %>% 
  knitr::kable(digits = 2)
season n FAC2 MB MGE NMB NMGE RMSE r COE IOA
2 spring (MAM) 1413 0.89 -5.64 17.15 -0.11 0.33 28.76 0.77 0.43 0.71
3 summer (JJA) 1185 0.83 3.08 9.20 0.13 0.38 12.16 0.38 0.03 0.51
1 autumn (SON) 1276 0.87 -0.74 15.37 -0.02 0.35 21.81 0.62 0.28 0.64
4 winter (DJF) 1945 0.80 7.07 21.68 0.13 0.40 31.64 0.71 0.34 0.67 
# 13 zmiennych
modStats(mydata = dane4_rf_test, 
         mod = "mod", 
         obs = "obs", 
         type = "season") %>% 
  knitr::kable(digits = 2)
season n FAC2 MB MGE NMB NMGE RMSE r COE IOA
2 spring (MAM) 1413 0.89 -5.82 16.94 -0.11 0.33 28.68 0.78 0.44 0.72
3 summer (JJA) 1185 0.85 2.04 8.74 0.08 0.36 11.58 0.42 0.08 0.54
1 autumn (SON) 1276 0.87 -0.93 15.35 -0.02 0.35 21.81 0.62 0.28 0.64
4 winter (DJF) 1945 0.78 9.23 22.86 0.17 0.42 32.24 0.71 0.30 0.65

Wnioski

Podsumowujac - wykonując to ćwiczenie poszerzyliśmy nasza wiedzę na temat przeprowadzania analizy danych metorologicznych, oraz o bazie danych GIOŚ. Poznaliśmy wiele możliwości skorzystania z tej wiedzy w celu przeprowadzenia jak najdokladniejszych analiz przy pomocy języka R. W przyszłości wykorystując tą wiedzę oraz zdobyte umiejętności będziemy w stanie osiagnąć zamierzone przez nas cele w stosunkowo krótkim czasie. Uważamy, że najważnieszym aspektem naszej pracy było zdobycie wiedzy z zakresu przeprowadzania prognozy stężenia PM10 przy pomocy RStudio.

Bibliografia

  1. A. Szulecka, R. Oleniacz, M. Rzeszutek (2017): Functionality of openair package in air pollution assessment and modeling – a case study of Krakow, Environmental Protection and Natural Resources, 28(2), 22-27. DOI: 10.1515/OSZN-2017-0009 (dostep 03.05.2021)
  2. J.N. Lott: The Quality Control of the Integrated Surface Hourly database (dostep 03.05.2021)
  3. D.C. Carslaw, K. Ropkins (2012): openair — An R package for air quality data analysis, Environmental Modelling & Software, 27–28(0), 52–61. DOI: 10.1016/j.envsoft.2011.09.008 (dostep 03.05.2021)
  4. D.B. Stephenson (2005): Data analysis methods in weather and climate research, on-line course (dostep 04.05.2021)
  5. https://bookdown.org/yihui/rmarkdown/html-document.html (dostep 05.05.2021)
  6. https://www.datadreaming.org/post/r-markdown-theme-gallery/ (dostep 05.05.2021)
  7. https://rpubs.com/danapower/577147 (dostep 05.05.2021)
  8. https://pl.wikipedia.org/wiki/R_(język_programowania) (dostep 10.05.2021)
  9. https://bookdown.org/yihui/rmarkdown/ (dostep 10.05.2021)
  10. https://dane.gov.pl/pl (dostep 15.05.2021)
  11. http://pbiecek.github.io/Przewodnik/Programowanie/jak_tworzyc_raporty.html (dostep 15.05.2021)
  12. https://tibble.tidyverse.org (dostep 15.05.2021)
  13. https://riptutorial.com/pl/r/example/2871/tworzenie-tabeli-danych (dostep 15.05.2021)
  14. https://cran.r-project.org/doc/contrib/wprowadzenie_do_R.pdf (dostep 15.05.2021)
  15. https://plotly.com/r/ (dostep 15.05.2021)
  16. https://cran.r-project.org/doc/contrib/Biecek-R-basics.pdf (dostep 15.05.2021)
  17. https://www.ncdc.noaa.gov/isd (dostep 15.05.2021)
  18. http://japoland.pl/blog/tokio-2/ (dostep 15.05.2021)
  19. https://powietrze.gios.gov.pl/pjp/content/about_us (dostep 16.05.2021)
  20. https://github.com/mlr-org/mlr/issues/1515 (dostęp 17.05.2021)
  21. https://community.rstudio.com/t/error-when-knitting-to-html-document-the-name-of-the-input-file-cannot-contain-the-special-shell-characters/45066/3 (dostep 17.05.2021)